For this assignment, we are going to investigate fires requiring the fire department to respond. Using data about the locations of firehouses and fires occurring in New York City, we want to know whether response times to fires differ across the city. Second, we will try to focus on one possible variable that could affect response times – the distance from the firehouse – and see whether we find the (expected) effect.
To keep this homework manageable, I am leaving out another part of the investigation: What is the effect of demographic and/or income characteristics of the neighborhood on response times. This is likely a bit more sensitive but also relevant from a public policy perspective.
We rely on two data sets.
NYC Open Data has data on all incidents responded to by fire companies. I have included the variable description file in the exercise folder. The following variables are available:
This dataset is only updated annually, and thus far only data from 2013 to 2018 is contained. The full dataset is also somewhat too large for an exercise (2.5M rows), so I suggest to limit yourself to a subset. I have added a file containing the subset of of only building fires (INCIDENT_TYPE_DESC == "111 - Building fire") for 2013 to 2018 only which yields about 14,000 incidents.
Unfortunately, the addresses of the incidents were not geocoded yet. Ideally, I would like you to know how to do this but am mindful about the hour or so required to get this done. So, here is the code. The geocodes (as far as they were returned successfully) are part of the data (as variables lat and lon).
library(ggmap)
Loading required package: ggplot2
Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df
Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
Please cite ggmap if you use it! See citation("ggmap") for details.
library(tidyverse)
[30m── [1mAttaching packages[22m ──────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✓[30m [34mtibble [30m 2.1.3 [32m✓[30m [34mpurrr [30m 0.3.3
[32m✓[30m [34mtidyr [30m 1.0.2 [32m✓[30m [34mdplyr [30m 0.8.3
[32m✓[30m [34mreadr [30m 1.3.1 [32m✓[30m [34mforcats[30m 0.4.0[39m
[30m── [1mConflicts[22m ─────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(dplyr)
# Open "building_fires" file
fire_building <- read_csv("data/building_fires.csv")
Parsed with column specification:
cols(
.default = col_character(),
IM_INCIDENT_KEY = [32mcol_double()[39m,
UNITS_ONSCENE = [32mcol_double()[39m,
TOTAL_INCIDENT_DURATION = [32mcol_double()[39m,
ZIP_CODE = [32mcol_double()[39m,
FIRE_ORIGIN_BELOW_GRADE_FLAG = [32mcol_double()[39m,
STORY_FIRE_ORIGIN_COUNT = [32mcol_double()[39m,
STANDPIPE_SYS_PRESENT_FLAG = [32mcol_double()[39m,
lon = [32mcol_double()[39m,
lat = [32mcol_double()[39m
)
See spec(...) for full column specifications.
NYC Open Data also provides data on the location of all 218 firehouses in NYC. Relevant for our analysis are the following variables: FacilityName, Borough, Latitude, Longitude
firehouses <- read_csv("data/FDNY_Firehouse_Listing.csv") %>%
dplyr::filter(!is.na(Latitude))
Parsed with column specification:
cols(
FacilityName = [31mcol_character()[39m,
FacilityAddress = [31mcol_character()[39m,
Borough = [31mcol_character()[39m,
Postcode = [32mcol_double()[39m,
Latitude = [32mcol_double()[39m,
Longitude = [32mcol_double()[39m,
`Community Board` = [32mcol_double()[39m,
`Community Council` = [32mcol_double()[39m,
`Census Tract` = [32mcol_double()[39m,
BIN = [32mcol_double()[39m,
BBL = [32mcol_double()[39m,
NTA = [31mcol_character()[39m
)
Note: 5 entries contain missing information, including on the spatial coordinates. We can exclude these for the exercise.
Provide a leaflet map of the highest severity fires (i.e. subset to the highest category in HIGHEST_LEVEL_DESC) contained in the file buiding_fires.csv. Ignore locations that fall outside the five boroughs of New York City. Provide at least three pieces of information on the incident in a popup.
# subset to the highest alarm
unique(fire_building$HIGHEST_LEVEL_DESC)
[1] "7 - Signal 7-5" "2 - 2nd alarm"
[3] "1 - More than initial alarm, less than Signal 7-5" "3 - 3rd alarm"
[5] "5 - 5th alarm" "4 - 4th alarm"
[7] NA "0 - Initial alarm"
[9] "75 - All Hands Working" "22 - Second Alarm"
[11] "55 - Fifth Alarm" "11 - First Alarm"
[13] "33 - Third Alarm" "44 - Fourth Alarm"
highest_alarm <- fire_building %>%
filter(HIGHEST_LEVEL_DESC=="7 - Signal 7-5"| HIGHEST_LEVEL_DESC=="75 - All Hands Working")
unique(highest_alarm$HIGHEST_LEVEL_DESC)
[1] "7 - Signal 7-5" "75 - All Hands Working"
head(highest_alarm)
library(leaflet)
Warning messages:
1: Unknown or uninitialised column: 'class'.
2: Unknown or uninitialised column: 'class'.
3: Unknown or uninitialised column: 'class'.
4: Unknown or uninitialised column: 'class'.
5: Unknown or uninitialised column: 'class'.
6: Unknown or uninitialised column: 'class'.
7: Unknown or uninitialised column: 'class'.
8: Unknown or uninitialised column: 'class'.
9: Unknown or uninitialised column: 'class'.
library(stringr)
highest_alarm_map <- leaflet(highest_alarm, options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles(provider = "Esri")%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 ) %>%
addCircleMarkers(lng = ~lon,lat = ~lat, radius = 1, popup = ~paste0("Date: ",str_extract(INCIDENT_DATE_TIME, pattern = "[0-9]+/[0-9]+/[0-9]+"), "<br/>", "Address: ", address, "<br/>", "Spread: ", str_sub(FIRE_SPREAD_DESC, 5,-1)))
highest_alarm_map
library(htmlwidgets)
saveWidget(highest_alarm_map, file="highest_alarm_map.html")
Start with the previous map. Now, distinguish the markers of the fire locations by PROPERTY_USE_DESC, i.e. what kind of property was affected. If there are too many categories, collapse some categories. Choose an appropriate coloring scheme to map the locations by type of affected property. Add a legend informing the user about the color scheme. Also make sure that the information about the type of affected property is now contained in the popup information. Show this map.
unique(highest_alarm$PROPERTY_USE_DESC)
[1] "579 - Motor vehicle or boat sales, services, repair"
[2] "429 - Multifamily dwelling"
[3] "419 - 1 or 2 family dwelling"
[4] "700 - Manufacturing, processing"
[5] "161 - Restaurant or cafeteria"
[6] "881 - Parking garage, (detached residential garage)"
[7] "564 - Laundry, dry cleaning"
[8] "400 - Residential, other"
[9] "559 - Recreational, hobby, home repair sales, pet store"
[10] "500 - Mercantile, business, other"
[11] "549 - Specialty shop"
[12] "331 - Hospital - medical or psychiatric"
[13] "131 - Church, mosque, synagogue, temple, chapel"
[14] "210 - Schools, non-adult, other"
[15] "519 - Food and beverage sales, grocery store"
[16] "449 - Hotel/motel, commercial"
[17] "173 - Bus station"
[18] "891 - Warehouse"
[19] "599 - Business office"
[20] "580 - General retail, other"
[21] "311 - 24-hour care Nursing homes, 4 or more persons"
[22] "UUU - Undetermined"
[23] "888 - Fire station"
[24] "900 - Outside or special property, other"
[25] "965 - Vehicle parking area"
[26] "100 - Assembly, other"
[27] "960 - Street, other"
[28] "332 - Hospices"
[29] "439 - Boarding/rooming house, residential hotels"
[30] "123 - Stadium, arena"
[31] "141 - Athletic/health club"
[32] "931 - Open land or field"
[33] "962 - Residential street, road or residential driveway"
[34] "899 - Residential or self-storage units"
[35] "322 - Alcohol or substance abuse recovery center"
[36] "170 - Passenger terminal, other"
[37] "557 - Personal service, including barber & beauty shops"
[38] "880 - Vehicle storage, other"
[39] "460 - Dormitory-type residence, other"
[40] "800 - Storage, other"
[41] "610 - Energy production plant, other"
[42] "926 - Outbuilding, protective shelter"
[43] "000 - Property Use, other"
[44] "511 - Convenience store"
[45] "162 - Bar or nightclub"
[46] "539 - Household goods, sales, repairs"
[47] "981 - Construction site"
[48] "180 - Studio/theater, other"
[49] "183 - Movie theater"
[50] "459 - Residential board and care"
[51] "241 - Adult education center, college classroom"
[52] "152 - Museum"
[53] "593 - Office: veterinary or research"
[54] "130 - Places of worship, funeral parlors, other"
[55] "160 - Eating, drinking places, other"
[56] "300 - Health care, detention, & correction, other"
[57] "642 - Electrical distribution"
[58] "963 - Street or road in commercial area"
[59] "342 - Doctor, dentist or oral surgeon office"
[60] "529 - Textile, wearing apparel sales"
[61] "321 - Mental retardation/development disability facility"
[62] "340 - Clinics, doctors offices, hemodialysis cntr, other"
[63] "882 - Parking garage, general vehicle"
[64] "150 - Public or government, other"
[65] "174 - Rapid transit station"
[66] "648 - Sanitation utility"
[67] "200 - Educational, other"
[68] "596 - Post office or mailing firms"
[69] "215 - High school/junior high school/middle school"
[70] "110 - Fixed-use recreation places, other"
[71] "592 - Bank"
[72] "365 - Police station"
[73] "571 - Service station, gas station"
[74] "182 - Auditorium, concert hall"
[75] "581 - Department or discount store"
[76] "808 - Outbuilding or shed"
[77] "839 - Refrigerated storage"
[78] "121 - Ballroom, gymnasium"
[79] "112 - Billiard center, pool hall"
[80] "NNN - None"
[81] "142 - Clubhouse"
[82] "140 - Clubs, other"
[83] "569 - Professional supplies, services"
[84] "124 - Playground"
[85] "363 - Reformatory, juvenile detention center"
[86] "974 - Aircraft loading area"
[87] "464 - Barracks, dormitory"
[88] "635 - Computer center"
[89] "211 - Preschool"
[90] "186 - Film/movie production studio"
[91] "181 - Live performance theater"
[92] "134 - Funeral parlor"
[93] "984 - Industrial plant yard - area"
[94] "629 - Laboratory or science lababoratory"
[95] "144 - Casino, gambling clubs"
[96] "936 - Vacant lot"
[97] "254 - Day care, in commercial property"
[98] "155 - Courthouse"
[99] "807 - Outside material storage area"
[100] "250 - Day care, other (Conversion only)"
[101] "615 - Electric-generating plant"
[102] "213 - Elementary school, including kindergarten"
[103] "143 - Yacht Club"
[104] "341 - Clinic, clinic-type infirmary"
[105] "639 - Communications center"
[106] "898 - Dock, marina, pier, wharf"
[107] "952 - Railroad yard"
There were 50 or more warnings (use warnings() to see the first 50)
highest_alarm %>%
group_by(PROPERTY_USE_DESC) %>%
count()%>%
arrange(desc (n))
library(stringr)
highest_alarm['class']=list(str_sub(highest_alarm$PROPERTY_USE_DESC,1,1))
highest_alarm$class[highest_alarm$class == "1"] <- "Assembly"
highest_alarm$class[highest_alarm$class == "2"] <- "Educational"
highest_alarm$class[highest_alarm$class == "3"] <- "Healthcare, Detention and Correction"
highest_alarm$class[highest_alarm$class == "4"] <- "Residential"
highest_alarm$class[highest_alarm$class == "5"] <- "Mercantile and Business"
highest_alarm$class[highest_alarm$class == "6"] <- "Energy Production Plant"
highest_alarm$class[highest_alarm$class == "7"] <- "Manufacturing and Processing"
highest_alarm$class[highest_alarm$class == "8"] <- "Storage"
highest_alarm$class[highest_alarm$class %in% c("9","U","N","0")] <- "Other Property"
highest_alarm$class <- factor(highest_alarm$class, levels= c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"))
highest_alarm %>%
group_by(class) %>%
count() %>%
arrange(desc(n))
library(RColorBrewer)
pal <- colorFactor(palette = "Set3", levels = c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"))
property_map <- leaflet(highest_alarm, options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles(provider = "CartoDB")%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 ) %>%
addCircleMarkers(lng = ~lon,lat = ~lat, radius = 1, color = ~pal(class), popup = ~paste0("<b>","Property Type: ", class,"</b>","<br/>","Date: ",str_extract(INCIDENT_DATE_TIME, pattern = "[0-9]+/[0-9]+/[0-9]+"), "<br/>", "Address: ", address, "<br/>", "Spread: ", str_sub(FIRE_SPREAD_DESC, 5,-1)))%>%
addLegend(pal = pal, values = c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"), opacity = 0.8, title = "Property Affected",position = "topleft")
property_map
saveWidget(property_map, file="property_map.html")
Add marker clustering, so that zooming in will reveal the individual locations but the zoomed out map only shows the clusters. Show the map with clusters.
cluster_property_map <- leaflet(highest_alarm, options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles(provider = "CartoDB")%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 ) %>%
addCircleMarkers(lng = ~lon,lat = ~lat, radius = 1, color = ~pal(class), clusterOptions = markerClusterOptions(), popup = ~paste0("<b>","Property Type: ", class,"</b>","<br/>","Date: ",str_extract(INCIDENT_DATE_TIME, pattern = "[0-9]+/[0-9]+/[0-9]+"), "<br/>", "Address: ", address, "<br/>", "Spread: ", str_sub(FIRE_SPREAD_DESC, 5,-1)))%>%
addLegend(pal = pal, values = c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"), opacity = 0.8, title = "Property Affected",position = "topleft")
cluster_property_map
saveWidget(cluster_property_map, file="cluster_property_map.html")
The second data file contains the locations of the 218 firehouses in New York City. Start with the non-clustered map (2b) and now adjust the size of the circle markers by severity (TOTAL_INCIDENT_DURATION or UNITS_ONSCENE seem plausible options). More severe incidents should have larger circles on the map. On the map, also add the locations of the fire houses. Add two layers (“Incidents”, “Firehouses”) that allow the user to select which information to show.
#highest_alarm$TOTAL_INCIDENT_DURATION)
summary(highest_alarm$TOTAL_INCIDENT_DURATION)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
870 4907 6803 7308 8918 428335 1
incidents_firehouse_map <-leaflet(highest_alarm, options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles(provider = "CartoDB")%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 ) %>%
addCircleMarkers(lng = ~lon,lat = ~lat, radius = ~TOTAL_INCIDENT_DURATION/5000 , color = ~pal(class), popup = ~paste0("<b>","Property Type: ", class,"</b>","<br/>","Date: ",str_extract(INCIDENT_DATE_TIME, pattern = "[0-9]+/[0-9]+/[0-9]+"), "<br/>", "Address: ", address, "<br/>", "Spread: ", str_sub(FIRE_SPREAD_DESC, 5,-1)), group = "Incidents")%>%
addLegend(pal = pal, values = c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"), opacity = 0.8, title = "Property Affected",position = "topleft")%>%
addMarkers(data = firehouses, lng = ~Longitude, lat = ~Latitude, group = "Firehouses")%>%
addLayersControl(overlayGroups = c("Incidents", "Firehouses"))
incidents_firehouse_map
saveWidget(incidents_firehouse_map, file="incidents_firehouse_map.html")
We now want to investigate whether the distance of the incident from the nearest firehouse varies across the city.
For all incident locations (independent of severity), identify the nearest firehouse and calculate the distance between the firehouse and the incident location. Provide a scatter plot showing the time until the first engine arrived (the variables INCIDENT_DATE_TIME and ARRIVAL_DATE_TIME) will be helpful.
library(sp)
library(sf)
Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(raster)
Attaching package: ‘raster’
The following object is masked from ‘package:dplyr’:
select
The following object is masked from ‘package:tidyr’:
extract
library(rgeos)
rgeos version: 0.5-2, (SVN revision 621)
GEOS runtime version: 3.7.2-CAPI-1.11.2
Linking to sp version: 1.3-1
Polygon checking: TRUE
# transform firehouse into sf object
firehouse_sf <- st_as_sf(firehouses, coords = c("Longitude", "Latitude"), crs = 4326)
# build an empty list to store the nearest distance
nearest <-list()
# iterate the fire_building and transform them into sp objects
# use st_distance to calculate the distances and find the nearest firehouse
for(i in 1:nrow(fire_building)){
point_sf <- st_as_sf(fire_building[i,], coords = c("lon", "lat"), crs = 4326)
nearest[i] <- min(st_distance(point_sf, firehouse_sf))
}
length(nearest)
[1] 14200
any(is.na(nearest))
[1] FALSE
# create a new column in fire_building about the nearest distance from firehouse
fire_building['distance'] <- unlist(nearest)
distance_and_time <- fire_building %>%
dplyr::select(INCIDENT_DATE_TIME, ARRIVAL_DATE_TIME, PROPERTY_USE_DESC,HIGHEST_LEVEL_DESC,BOROUGH_DESC, lon, lat, distance)
head(distance_and_time)
incident_time <- as.POSIXct(strptime(distance_and_time[['INCIDENT_DATE_TIME']], format = "%m/%d/%Y %H:%M:%S %p"))
arrival_time <- as.POSIXct(strptime(distance_and_time[['ARRIVAL_DATE_TIME']], format = "%m/%d/%Y %H:%M:%S %p"))
waiting_time <- arrival_time-incident_time
distance_and_time['waiting_time_secs'] <- as.numeric(waiting_time)
#There are some mistakes in original records like wrong AM/PM
#Transfrom those wrong records (negative numbers) by adding 12 hours back
for(i in 1:nrow(distance_and_time)){
if (!is.na(distance_and_time[i,'waiting_time_secs'])&distance_and_time[i,'waiting_time_secs'] < 0) {
distance_and_time[i,'waiting_time_secs'] <- distance_and_time[i,'waiting_time_secs']+12*60*60
}
}
head(distance_and_time)
summary(distance_and_time)
INCIDENT_DATE_TIME ARRIVAL_DATE_TIME PROPERTY_USE_DESC HIGHEST_LEVEL_DESC BOROUGH_DESC lon
Length:14200 Length:14200 Length:14200 Length:14200 Length:14200 Min. :-74.25
Class :character Class :character Class :character Class :character Class :character 1st Qu.:-73.97
Mode :character Mode :character Mode :character Mode :character Mode :character Median :-73.92
Mean :-73.92
3rd Qu.:-73.87
Max. :-73.09
lat distance waiting_time_secs
Min. :40.50 Min. : 7.75 Min. : 12.0
1st Qu.:40.67 1st Qu.: 362.06 1st Qu.: 172.0
Median :40.72 Median : 561.04 Median : 208.0
Mean :40.73 Mean : 666.23 Mean : 421.1
3rd Qu.:40.80 3rd Qu.: 863.99 3rd Qu.: 250.0
Max. :41.60 Max. :101812.17 Max. :86816.0
NA's :23
library(ggplot2)
library(ggthemes)
library(plotly)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Attaching package: ‘plotly’
The following object is masked from ‘package:raster’:
select
The following object is masked from ‘package:ggmap’:
wind
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
df_plot <-distance_and_time %>%
filter(!is.na(waiting_time_secs))%>%
filter(distance < 5000) %>%
filter(waiting_time_secs < 5000)
p <-ggplot(df_plot, aes(x = distance, y = waiting_time_secs))+
geom_point(alpha = 0.5)+
labs(x = "Distance From The Nearest Firehouse (m)", y = "Waiting Time For The First Engine (secs)")+
theme_clean()
p
Now also visualize the patterns separately for severe and non-severe incidents (use HIGHEST_LEVEL_DESC but feel free to reduce the number of categories). What do you find?
fire_levels <- distance_and_time%>%
filter(!is.na(HIGHEST_LEVEL_DESC))
fire_levels['level']=list(str_sub(fire_levels$HIGHEST_LEVEL_DESC,1,2))
fire_levels$level[fire_levels$level %in% c("7 ", "75")] <- "High Alarm"
fire_levels$level[fire_levels$level %in% c("5 ", "55", "4 ", "44", "3 ", "33")] <- "Medium Alarm"
fire_levels$level[fire_levels$level %in% c("2 ", "22", "11", "0 ")] <- "Low Alarm"
fire_levels$level[fire_levels$level %in% c("1 ")] <- "Undefined Alarm"
fire_levels$level <- factor(fire_levels$level, levels= c("Low Alarm","Medium Alarm","High Alarm","Undefined Alarm"))
fire_levels_plot <- fire_levels %>%
filter(!is.na(waiting_time_secs))%>%
filter(distance < 5000) %>%
filter(waiting_time_secs < 1000)
ggplot(fire_levels_plot, aes(x=waiting_time_secs, y=distance, color = level))+
geom_point(alpha=0.4)+
facet_grid(~ level)+
labs(x = "Waiting Time For The First Engine (secs)", y = "Distance From The Nearest Firehouse (m)")+
theme_bw()
Provide a map visualization of response times. Investigate whether the type of property affected (PROPERTY_USE_DESC) or fire severity (HIGHEST_LEVEL_DESC) play a role here.
response_time <- fire_levels
response_time['class']=list(str_sub(response_time$PROPERTY_USE_DESC,1,1))
response_time$class[response_time$class == "1"] <- "Assembly"
response_time$class[response_time$class == "2"] <- "Educational"
response_time$class[response_time$class == "3"] <- "Healthcare, Detention and Correction"
response_time$class[response_time$class == "4"] <- "Residential"
response_time$class[response_time$class == "5"] <- "Mercantile and Business"
response_time$class[response_time$class == "6"] <- "Energy Production Plant"
response_time$class[response_time$class == "7"] <- "Manufacturing and Processing"
response_time$class[response_time$class == "8"] <- "Storage"
response_time$class[response_time$class %in% c("9","U","N","0")] <- "Other Property"
response_time$class <- factor(response_time$class, levels= c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"))
head(response_time)
low <- filter(response_time, level == "Low Alarm")
medium <- filter(response_time, level == "Medium Alarm")
high <- filter(response_time, level == "High Alarm")
pal_level <- colorFactor(palette = c("#f03b20", "#feb24c", "#ffeda0"),
levels = c("High Alarm","Medium Alarm","Low Alarm"))
alarm_map <- leaflet(options = leafletOptions(minZoom = 8, dragging = TRUE)) %>%
addProviderTiles("CartoDB.DarkMatter",options = providerTileOptions(attribution = ""))%>%
addCircleMarkers(data = low, lng = ~lon, lat = ~lat, radius = ~waiting_time_secs/5000, fillOpacity=0.3,
color = ~pal_level(level), group = "Low Alarm", popup=~paste("Alarm Level: ",level,
"<br>Response Time: ", waiting_time_secs," seconds", "<br>Property Type: ", class)) %>%
addCircleMarkers(data = medium, lng = ~lon, lat = ~lat, radius = ~waiting_time_secs/5000, fillOpacity=0.3,
color = ~pal_level(level), group = "Medium Alarm", popup=~paste("Alarm Level: ",level,
"<br>Response Time: ", waiting_time_secs, " seconds", "<br>Property Type: ", class)) %>%
addCircleMarkers(data = high, lng = ~lon, lat = ~lat, radius = ~waiting_time_secs/5000, fillOpacity=0.3,
color = ~pal_level(level), group = "High Alarm", popup=~paste("Alarm Level: ",level,
"<br>Response Time: ", waiting_time_secs, " seconds", "<br>Property Type: ", class)) %>%
setView(lat= 40.712742, lng=-74.013382, zoom = 10) %>%
addLegend(pal = pal_level, values = c("High Alarm","Medium Alarm", "Low Alarm"), opacity = 0.7, title = "Alarm Level",
position = "topleft")%>%
addLayersControl(overlayGroups = c("High Alarm","Medium Alarm","Low Alarm"))
alarm_map
saveWidget(alarm_map, file="alarm_map.html")
long_response <- subset(response_time, waiting_time_secs > 500)
fireIcons <- icons(
iconUrl = "data/redflame.png",
iconWidth = 15, iconHeight = 15,
iconAnchorX = 7.5, iconAnchorY = 8.5
)
pal_class <- colorFactor(palette = "Tableau10", levels = c("Assembly","Educational","Healthcare, Detention and Correction", "Residential", "Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"))
response_residential <- leaflet(options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles(provider = "CartoDB") %>%
addCircleMarkers(data = response_time, lng = ~lon, lat = ~lat, radius = ~waiting_time_secs/5000, color = ~pal_class(class),
fillOpacity=0.3, popup = ~paste("Property Type: ", class, "<br>Response Time: ", waiting_time_secs," seconds"))%>%
addMarkers(data= long_response,icon = fireIcons,
popup = ~paste("Property Type: ", class, "<br>Response Time: ", waiting_time_secs," seconds"))%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 )
Assuming "lon" and "lat" are longitude and latitude, respectively
response_residential
#saveWidget(response_residential, file="response_residential.html")
According to the plot above, those big circles and fire flame icons are the incidents with a long response time.Clicking the popup, we can see most of them are residential properties. Then we remove those “Residential” records and find some incidents happened in businesses and assembly also had a long response time.
fireIcons2 <- icons(
iconUrl = "data/flame.png",
iconWidth = 15, iconHeight = 15,
iconAnchorX = 7.5, iconAnchorY = 8.5
)
pal_property <- colorFactor(palette = "Spectral", levels = c("Assembly","Educational","Healthcare, Detention and Correction","Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"))
response_property_map <- leaflet(options = leafletOptions(minZoom = 5, dragging = TRUE))%>%
addProviderTiles("Esri.WorldImagery", options = providerTileOptions(attribution = "")) %>%
addCircleMarkers(data = subset(response_time, response_time$class != "Residential"), lng = ~lon, lat = ~lat, radius = ~waiting_time_secs/5000, color = ~pal_property(class), fillOpacity=0.7, popup = ~paste("Property Type: ", class, "<br>Response Time: ", waiting_time_secs," seconds"))%>%
addMarkers(data= subset(long_response, long_response$class != "Residential"),icon = fireIcons2,
popup = ~paste("Property Type: ", class, "<br>Response Time: ", waiting_time_secs," seconds"))%>%
addLegend(pal = pal_property, values = c("Assembly","Educational","Healthcare, Detention and Correction","Mercantile and Business", "Energy Production Plant", "Manufacturing and Processing", "Storage", "Other Property"), opacity = 0.7, title = "Property Affected",position = "topleft")%>%
setView( lat= 40.712742, lng=-74.013382, zoom = 10 )
Assuming "lon" and "lat" are longitude and latitude, respectively
response_property_map
#saveWidget(response_property_map, file="response_property_map.html")
Show a faceted choropleth map indicating how response times have developed over the years. What do you find?
response_borough<-distance_and_time %>%
select(INCIDENT_DATE_TIME,BOROUGH_DESC, waiting_time_secs)
response_borough['borough'] = list(str_sub(response_borough$BOROUGH_DESC,1,1))
response_borough$borough[response_borough$borough == "1"] <- "Manhattan"
response_borough$borough[response_borough$borough == "2"] <- "Bronx"
response_borough$borough[response_borough$borough == "3"] <- "Staten Island"
response_borough$borough[response_borough$borough == "4"] <- "Brooklyn"
response_borough$borough[response_borough$borough == "5"] <- "Queens"
response_borough$borough <- factor(response_borough$borough, levels= c("Manhattan","Bronx","Staten Island","Brooklyn","Queens"))
response_borough['year'] = list(str_sub(response_borough$INCIDENT_DATE_TIME,7,10))
response_borough$year <- factor(response_borough$year, levels= c("2013","2014","2015","2016","2017","2018"))
head(response_borough)
subset(response_borough, response_borough$borough == "Queens"&response_borough$year =="2013")
average_response_time <-response_borough %>%
filter(!is.na(waiting_time_secs))%>%
group_by(borough, year) %>%
summarise(mean_response_time = round(mean(waiting_time_secs),2))
average_response_time
library(rgdal)
rgdal: version: 1.4-8, (SVN revision 845)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE
Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
Linking to sp version: 1.3-2
borough <- readOGR("data/borough_boundaries.geojson", verbose=FALSE)
borough@data
shp_response <- borough@data %>%
right_join(average_response_time, by = c("boro_name"= "borough"))
Column `boro_name`/`borough` joining factors with different levels, coercing to character vector
shp_response
borough@data <-shp_response %>%
filter(year =="2013")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
231.6 351.8 387.8 393.4 497.8 498.3
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2013 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2013",
position = "topleft", opacity=0.7)
borough@data <-shp_response %>%
filter(year =="2014")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
234.8 297.1 380.4 465.0 492.6 920.3
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2014 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2014",
position = "topleft", opacity=0.7)
borough@data <-shp_response %>%
filter(year =="2015")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
188.1 217.4 396.9 336.6 398.4 482.6
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2015 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2013",
position = "topleft", opacity=0.7)
borough@data <-shp_response %>%
filter(year =="2016")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
226.1 233.4 359.7 353.4 410.7 537.0
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2016 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2016",
position = "topleft", opacity=0.7)
borough@data <-shp_response %>%
filter(year =="2017")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
223.3 228.8 427.3 410.8 535.7 638.9
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2017 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2017",
position = "topleft", opacity=0.7)
borough@data <-shp_response %>%
filter(year =="2018")
summary(borough$mean_response_time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
236.0 244.1 246.9 396.6 248.8 1007.0
pal_response <- colorNumeric("PuOr", domain = borough$mean_response_time)
map2018 <-borough %>%
leaflet()%>%
addProviderTiles("CartoDB")%>%
addPolygons(weight = 1, color = ~pal_response(mean_response_time), fillOpacity = 1,
label = ~paste0("Mean Response Time: ", mean_response_time, "seconds"),
highlightOptions = highlightOptions(weight = 5, color = "white", bringToFront = TRUE))%>%
addLegend(pal = pal_response, values = ~ mean_response_time, title = "2018",
position = "topleft", opacity=0.7)
library(mapview)
facet_map <- sync(map2013, map2014, map2015, map2016, map2017,map2018, ncol = 3, sync = "all")
'mapview::sync' is deprecated.
Use 'leafsync::sync' instead.
See help("Deprecated") and help("leafsync-deprecated").'mapview::latticeView' is deprecated.
Use 'leafsync::latticeView' instead.
See help("Deprecated") and help("mapview-deprecated").
facet_map
Please follow the instructions to submit your homework. The homework is due on Wednesday, March 25.
If you do come across something online that provides part of the analysis / code etc., please no wholesale copying of other ideas. We are trying to evaluate your abilities to visualize data, not the ability to do internet searches. Also, this is an individually assigned exercise – please keep your solution to yourself.